/**
 * @file LinearInterpolation.h
 * @brief The implementation of Linear interpolation operator.
 * @author XDDDD
 * @version 
 * @date 2021-06-10
 */

#ifndef __PAISLEYPARK__LINEARINTERPOLATION_H__
#define __PAISLEYPARK__LINEARINTERPOLATION_H__

#include "InterpolationOperator.h"

template <size_t DIM>
class LinearInterpolation : public InterpolationOperator<DIM> {
public:
	void operator()(MultigridVector<DIM>& _v, const size_t _M) override;
};

template <size_t DIM>
void LinearInterpolation<DIM>::operator()(MultigridVector<DIM>& _v, const size_t _M) {
	if(_M == 0 || _M >= _v.get_N()) {
		std::cout << "Invalid _M" << std::endl;
		std::exit(-1);
	}
	std::vector<size_t> coord_2h;
	std::vector<size_t> coord_h;
	coord_2h.clear();
	coord_2h.resize(DIM);
	coord_h.clear();
	coord_h.resize(DIM);
	int d[DIM] = {0};
	size_t pole = 0;
	size_t num_1d = (1 << _M) + 1;
	double coef = 1.0;
	bool is_bound = false;
	//clear.
	_v.cleanpart(_M + 1);
	//Go through all the nodes in the coareser grid.
	while(pole != DIM) {
		while(pole != DIM) {
			pole = 0;
			coef = 1.0;
			is_bound = false;
			for(size_t i = 0; i != DIM; i++) {
				if(d[i] != 0)
					coef *= 0.5;
				coord_h[i] = d[i] + 2*coord_2h[i];		//May exist bug.
			}
			for(size_t i = 0; i != DIM; i++) {
				if(!(coord_h[i] > 0 && coord_h[i] < 2*num_1d - 2)) {
					is_bound = true;
				}
			}
			if(!is_bound) {	
				_v(coord_h, _M + 1) += coef*_v(coord_2h, _M);
			}
			d[pole] = -2*d[pole] - 1;
			while(d[pole] == -3) {
				d[pole] = 0;
				pole++;
				if(pole == DIM)
					break;
				d[pole] = -2*d[pole] - 1;
			}
		}
		pole = 0;
		coord_2h[pole]++;
		while(coord_2h[pole] == num_1d) {
			coord_2h[pole] = 0;
			pole++;
			if(pole == DIM)
				break;
			coord_2h[pole]++;
		}
	}
};

#else
//Do nothing.
#endif
